function Y = compute_D(X, n, p)
	Y = zeros(n, n*(p+1));
	for i=0:p
		Y(1:n,1:n) = Y(1:n,1:n) + X((1:n)+n*i,(1:n)+n*i);
	end
	for k=1:p
		for i=0:(p-k)
			Y(1:n,(1:n)+n*k) = Y(1:n,(1:n)+n*k) + 2*X((1:n)+n*i,(1:n)+n*(i+k));
		end
	end
end
